Who Benefits from Religious Attendance?

Heterogeneous Causal Effects on Well-being and Cooperation in New Zealand

Joseph A. Bulbulia
joseph.bulbulia@vuw.ac.nz

Victoria University of Wellington, New Zealand

Causal Effects Effects of Religious Service

Three Questions:

  1. Flourishing – does attendance enhance multi‑dimensional well‑being (“the good life”)?

  2. Cooperation – does attendance spur pro-social behaviour?

  3. Individual Differences - who responds differently?

Why Causality Matters

On Average The Walk Is Flat

Method

  1. Track same people over 6 years (n = 46,377)
  2. Compare those who change attendance vs those who don’t
  3. Use Machine Learning to find who benefits most
  4. Account for selection bias and dropout
  5. Focus on causal not correlational effects
  • Probability panel (\approx 1% of electorate; 70–80% annual retention).
  • Domains: personality, social attitudes, values, religion, employment, prejudice, religion \dots osf-link
  • Multidisciplinary team > 60 researchers
  • Responded at Wave 10 (2018/19) baseline
  • Reported service attendance at baseline (Yes/No)
  • Inverse‑probability of censoring weights adjust attrition
  • “Censored” if data collected during national lockdown (Wave 10)
  • Final analytic sample: 46,377 adults.

Note

Technical details in Supplements

Stage Tool Purpose
1. Target Causal estimand definition Binary & shift interventions
2. Learn SuperLearner (semi-parametric ML) Outcome & censoring models
3. Adjust IPC weights + census post-stratification Population inference
4. Validate 10-fold cross-validation Out-of-sample performance
5. Subgroup Causal forests & policy trees Detect heterogeneity

Note

Technical details in Supplements

Estimand Description Example
Binary All or nothing 0 vs 1 service/month
Shift-up Increase unless weekly 0 → 1, 1 → 2, 2 → 3, 3 → 4
Shift-down Decrease to zero 4 → 0, 3 → 0, 2 → 0, 1 → 0

Note

Technical details in Supplements

Daily Survey Throughput

Study 1 | Flourishing

1.1 Binary intervention (3 waves)

Figure 1: ATE: Three-Wave: Binary Intervention

1.2 Shift‑up vs binary (3 waves)

Figure 2: ATE: Three-wave Shift Intervention

1.3 Sensitivity: socialising hours (3 waves)

Figure 3: ATE: Three-wave Soft Intervention

1.4 Duration: 3 vs 6 waves

Figure 4: ATE: Six-wave Soft Intervention GAIN: +1

1.5 Gain vs loss (6 waves)

Figure 5: ATE: Six-wave Soft Intervention LOSS: -1

1.6 Good Life: Compare 6 waves WEEKLY/NULL vs 6 waves ZERO/NULL

Figure 6: CATE RATE

1.7 Findings – Flourishing

  • Meaning & purpose show the most robust positive effect.
  • Binary gains detectable; shift‑based gains weaker but present.
  • Hour‑for‑hour socialising does not explain the effect.
  • Sexual satisfaction follows a similar pattern (consistent with prior literature)
Figure 7: CATE rate

Study 2 | Cooperation

2.1 Binary intervention (3 waves)

Figure 8: ATE: Three-Wave: Binary Intervention

2.2 Shift‑up vs binary (3 waves)

Figure 9: ATE: Three-wave Shift Intervention

2.3 Gain vs loss (3 waves)

Figure 10: ATE: Three-wave Soft Intervention

2.4 Duration: 3 vs 6 waves (gain)

Figure 11: ATE: Six-wave Soft Intervention WEEKLY/ZER0

2.5 Gain vs loss (6 waves)

Figure 12: ATE: Six-wave Soft Intervention LOSS: -1

2.6 Findings – Cooperation

  • Gains in attendance cause small but reliable boosts to donations & neighbourliness.
  • Losses erode cooperation more than gains enhance it (asymmetry).
  • Perceived social support and belonging are not clearly affected
Figure 13: CATE RATE

Study 3a | Personalised Flourishing

Forgiveness

Figure 14: Policy Trsees: Forgiveness

Gratitude

Figure 15: Policy Trees: Gratitude

Fatigue

Figure 16: Policy Trees for Hlth Fatigue Outcome

Meaning Sense

Figure 17: Policy Trees: Meaning Sense

Sense of Neighbourhood Community

Figure 18: Policy Trees: Neighbourhood Community

Self Esteem

Figure 19: Policy Trees: Self Esteem

Study 3b Personalised Cooperation

Donations

Figure 20: Policy Trees for log Charity Donate Outcome

3. Policy Trees Summary

Effect-Modifiers that Matter

  • Age, Fatigue, Sleep, Big-6 personality, Perfectionism, Work

  • Economic strain and hours of unpaid housework

  • Not major moderators: gender (rather housework/income) or ethnicity (rather personality)

Summary Graphs

Figure 21: Policy Trsees: Forgiveness
Figure 22: Policy Trees: Gratitude
Figure 23: Policy Trees for Hlth Fatigue Outcome
Figure 24: Policy Trees: Meaning Sense
Figure 25: Policy Trees: Neighbourhood Community
Figure 26: Policy Trees: Self Esteem
Figure 27: Policy Trees for log Charity Donate Outcome

Summary

  1. Signals modest + noisy (NZ context).
  2. Attendance mainly ↑ meaning ↑ prosociality.
  3. Effect-modification: virtues/personality/demographics, not “identities” (gender/ethnicity/etc.)

Thanks

  • 76,409 individuals who have participated in the New Zealand Attitudes and Values Study since 2009.
  • Templeton Religion Trust
  • University of Auckland
  • Victoria University
  • Georgia State University
  • Don E. Davis
  • Ken Rice
  • Geoff Troughton
  • Chris G. Sibley
  • Grad Students in the EPIC Lap
  • & other collaborators

Thank You!

S1 Daily Data Collection

Figure 28 presents the New Zealand Attitudes and Values Study Data Collection (2018 retained cohort) from years 2018-2014 (NZAVS time 10–time 15).

Figure 28: Historgram of New Zealand Attitudes and Values Study Daily Data Collection for Time 10 cohort: years 2018-2024.

S2: Measures and Demographic Statistics

Measures

Baseline Covariates

Age

We asked participants’ ages in an open-ended question (“What is your age?” or “What is your date of birth”) (Chris G. Sibley 2021).

Items:

  • What is your date of birth?

Agreeableness

Mini-IPIP6 Agreeableness dimension: (i) I sympathize with others’ feelings. (ii) I am not interested in other people’s problems. (r) (iii) I feel others’ emotions. (iv) I am not really interested in others. (r) (Chris G. Sibley et al. 2011).

Items:

  • I sympathize with others’ feelings.
  • I am not interested in other people’s problems.
  • I feel others’ emotions.
  • I am not really interested in others (reversed).

Alcohol Frequency

Participants could chose between the following responses: ‘(1 = Never - I don’t drink, 2 = Monthly or less, 3 = Up to 4 times a month, 4 = Up to 3 times a week, 5 = 4 or more times a week, 6 = Don’t know)’ (Health 2013).

Items:

  • “How often do you have a drink containing alcohol?”

Alcohol Intensity

Participants responded using an open-ended box (Health 2013).

Items:

  • “How many drinks containing alcohol do you have on a typical day when drinking alcohol? (number of drinks on a typical day when drinking)”

Social Belonging

We assessed felt belongingness with three items adapted from the Sense of Belonging Instrument (Hagerty & Patusky, 1995): (1) “Know that people in my life accept and value me”; (2) “Feel like an outsider”; (3) “Know that people around me share my attitudes and beliefs”. Participants responded on a scale from 1 (Very Inaccurate) to 7 (Very Accurate). The second item was reversely coded (Hagerty and Patusky 1995).

Items:

  • Know that people in my life accept and value me.
  • Feel like an outsider (reversed).
  • Know that people around me share my attitudes and beliefs.

Born in NZ

Coded binary (1 = New Zealand; 0 = elsewhere.) (Chris G. Sibley 2021).

Items:

  • Where were you born? (please be specific, e.g., which town/city?)

Conscientiousness

Mini-IPIP6 Conscientiousness dimension: (i) I get chores done right away. (ii) I like order. (iii) I make a mess of things. (r) (iv) I often forget to put things back in their proper place. (r) (Chris G. Sibley et al. 2011).

Items:

  • I get chores done right away.
  • I like order.
  • I make a mess of things.
  • I often forget to put things back in their proper place.

Education Level

We asked participants, “What is your highest level of qualification?”. We coded participans highest finished degree according to the New Zealand Qualifications Authority. Ordinal-Rank 0-10 NZREG codes (with overseas school qualifications coded as Level 3, and all other ancillary categories coded as missing) (Chris G. Sibley 2021).

Items:

  • What is your highest level of qualification?

Employed (Binary)

Binary response: (0 = No, 1 = Yes) (Statistics New Zealand 2017).

Items:

  • Are you currently employed (This includes self-employed of casual work)?

Ethnicity

Coded string: (1 = New Zealand European; 2 = Māori; 3 = Pacific; 4 = Asian) (Statistics New Zealand 2017).

Items:

  • Which ethnic group(s) do you belong to?

Extraversion

Mini-IPIP6 Extraversion dimension: (i) I am the life of the party. (ii) I don’t talk a lot. (r) (iii) I keep in the background. (r) (iv) I talk to a lot of different people at parties (Chris G. Sibley et al. 2011).

Items:

  • I am the life of the party.
  • I don’t talk a lot (reversed).
  • I keep in the background (reversed).
  • I talk to a lot of different people at parties.

Hlth Disability Binary

We assessed disability with a one-item indicator adapted from Verbrugge (1997). It asks, “Do you have a health condition or disability that limits you and that has lasted for 6+ months?” (1 = Yes, 0 = No) (Verbrugge 1997).

Items:

  • Do you have a health condition or disability that limits you and that has lasted for 6+ months?

Honesty Humility

Mini-IPIP6 Honesty-Humility dimension: (i) I feel entitled to more of everything. (r) (ii) I deserve more things in life. (r) (iii) I would like to be seen driving around in a very expensive car. (r) (iv) I would get a lot of pleasure from owning expensive luxury goods. (r) (Chris G. Sibley et al. 2011).

Items:

  • I feel entitled to more of everything (reversed).
  • I deserve more things in life (reversed).
  • I deserve more things in life (reversed).
  • I would get a lot of pleasure from owning expensive luxury goods (reversed).

log Hours Children

We took the natural log of the response + 1 (Chris G. Sibley et al. 2011).

Items:

  • Hours spent…looking after children.

log Hours Community

No information available for this variable.

log Hours Commute

We took the natural log of the response + 1 (Chris G. Sibley 2021).

Items:

  • Hours spent…travelling/commuting.

log Hours of Exercise

We took the natural log of the response + 1 (Chris G. Sibley et al. 2011).

Items:

  • Hours spent…exercising/physical activity.

log Hours Housework

We took the natural log of the response + 1 (Chris G. Sibley et al. 2011).

Items:

  • Hours spent…housework/cooking.

log Household Income

We took the natural log of the response + 1 (Chris G. Sibley 2021).

Items:

  • Please estimate your total household income (before tax) for the year XXXX.

Male (Binary)

Here, we coded all those who responded as Male as 1, and those who did not as 0 (Fraser et al. 2020).

Items:

  • We asked participants’ gender in an open-ended question: “what is your gender?”

Neuroticism

Mini-IPIP6 Neuroticism dimension: (i) I have frequent mood swings. (ii) I am relaxed most of the time. (r) (iii) I get upset easily. (iv) I seldom feel blue. (r) (Chris G. Sibley et al. 2011).

Items:

  • I have frequent mood swings.
  • I am relaxed most of the time (reversed).
  • I get upset easily.
  • I seldom feel blue (reversed).

Not Heterosexual Binary

Open-ended question, coded as binary (not heterosexual = 1) (Greaves et al. 2017).

Items:

  • How would you describe your sexual orientation? (e.g., heterosexual, homosexual, straight, gay, lesbian, bisexual, etc.)

NZ Deprevation Index 2018

Numerical: (1-10) (Atkinson, Salmond, and Crampton 2019).

Items:

  • New Zealand Deprivation - Decile Index - Using 2018 Census Data

NZSEI (Occupational Prestige Index)

This index uses the income, age, and education of a reference group, in this case, the 2013 New Zealand census, to calculate a score for each occupational group. Scores range from 10 (Lowest) to 90 (Highest). This list of index scores for occupational groups was used to assign each participant a NZSEI-13 score based on their occupation (Fahy, Lee, and Milne 2017).

Items:

  • We assessed occupational prestige and status using the New Zealand Socio-economic Index 13 (NZSEI-13).

Openness

Mini-IPIP6 Openness to Experience dimension: (i) I have a vivid imagination. (ii) I have difficulty understanding abstract ideas. (r) (iii) I do not have a good imagination. (r) (iv) I am not interested in abstract ideas. (r) (Chris G. Sibley et al. 2011).

Items:

  • I have a vivid imagination.
  • I have difficulty understanding abstract ideas (reversed).
  • I do not have a good imagination (reversed).
  • I am not interested in abstract ideas (reversed).

Parent (Binary)

Parents were coded as 1, while the others were coded as 0 (Chris G. Sibley 2021).

Items:

  • If you are a parent, in which year was your eldest child born?

Partner Binary

Coded as binary (has partner = 1) (Chris G. Sibley 2021).

Items:

  • What is your relationship status? (e.g., single, married, de-facto, civil union, widowed, living together, etc.)

Political Conservative

  • Please rate how politically liberal versus conservative you see yourself as being.

Rural Gch 2018 Levels

“Participants residence locations were coded according to a five-level ordinal categorisation ranging from Urban to Rural.” (Whitehead et al. 2023).

Items:

  • High Urban Accessibility = 1, Medium Urban Accessibility = 2, Low Urban Accessibility = 3, Remote = 4, Very Remote = 5.

Right Wing Authoritarianism

Right Wing Authoritarianism was measured using the following items (Altemeyer 1996).

Items:

  • It is always better to trust the judgment of the proper authorities in government and religion than to listen to the noisy rabble-rousers in our society who are trying to create doubt in people’s minds.
  • It would be best for everyone if the proper authorities censored magazines so that people could not get their hands on trashy and disgusting material.
  • Our country will be destroyed some day if we do not smash the perversions eating away at our moral fibre and traditional beliefs.
  • People should pay less attention to The Bible and other old traditional forms of religious guidance, and instead develop their own personal standards of what is moral and immoral.
  • Atheists and others who have rebelled against established religions are no doubt every bit as good and virtuous as those who attend church regularly.
  • Some of the best people in our country are those who are challenging our government, criticizing religion, and ignoring the “normal way” things are supposed to be done (reversed).

Sample Frame Opt in (Binary)

Code string (Binary): (0 = No, 1 = Yes) (Chris G. Sibley 2021).

Items:

  • Participant was not randomly sampled from the New Zealand Electoral Roll.

Social Dominance Orientation

Social Dominance Orientation was measured using the following items (Sidanius and Pratto 1999).

Items:

  • It is OK if some groups have more of a chance in life than others.
  • Inferior groups should stay in their place.
  • To get ahead in life, it is sometimes okay to step on other groups.
  • We should have increased social equality (reversed).
  • It would be good if groups could be equal (reversed).
  • We should do what we can to equalise conditions for different groups (reversed).

Smoker (Binary)

Binary smoking indicator (0 = No, 1 = Yes) (Chris G. Sibley 2021).

Items:

  • Do you currently smoke tobacco cigarettes?

Sample Demographic Statistics

Table 1 presents sample demographic statistics.

Table 1: Demographic statistics for New Zealand Attitudes and Values Study Cohort in the baseline weave (NZAVS time 10, spanning years 2018-2019).
Baseline Wave
(N=46377)
Age
Mean (SD) 48.6 (13.9)
Median [Min, Max] 51.0 [18.0, 99.0]
Agreeableness
Mean (SD) 5.35 (0.988)
Median [Min, Max] 5.50 [1.00, 7.00]
Missing 400 (0.9%)
Alcohol Frequency
Mean (SD) 2.16 (1.34)
Median [Min, Max] 2.00 [0, 5.00]
Missing 1598 (3.4%)
Alcohol Intensity
Mean (SD) 2.17 (2.17)
Median [Min, Max] 2.00 [0, 30.0]
Missing 2751 (5.9%)
Social Belonging
Mean (SD) 5.14 (1.08)
Median [Min, Max] 5.33 [1.00, 7.00]
Missing 397 (0.9%)
Born in NZ
0 10041 (21.7%)
1 36184 (78.0%)
Missing 152 (0.3%)
Conscientiousness
Mean (SD) 5.11 (1.06)
Median [Min, Max] 5.25 [1.00, 7.00]
Missing 392 (0.8%)
Education Level
no_qualification 1177 (2.5%)
cert_1_to_4 16277 (35.1%)
cert_5_to_6 5821 (12.6%)
university 12311 (26.5%)
post_grad 5034 (10.9%)
masters 3857 (8.3%)
doctorate 1111 (2.4%)
Missing 789 (1.7%)
Employed (binary)
0 9501 (20.5%)
1 36865 (79.5%)
Missing 11 (0.0%)
Ethnicity
euro 36915 (79.6%)
maori 5311 (11.5%)
pacific 1109 (2.4%)
asian 2453 (5.3%)
Missing 589 (1.3%)
Extraversion
Mean (SD) 3.91 (1.20)
Median [Min, Max] 4.00 [1.00, 7.00]
Missing 392 (0.8%)
Hlth Disability Binary
Mean (SD) 0.225 (0.418)
Median [Min, Max] 0 [0, 1.00]
Missing 893 (1.9%)
Honesty Humility
Mean (SD) 5.41 (1.18)
Median [Min, Max] 5.50 [1.00, 7.00]
Missing 396 (0.9%)
Hours Children
Mean (SD) 14.0 (32.3)
Median [Min, Max] 0 [0, 168]
Missing 1442 (3.1%)
Hours Commute
Mean (SD) 5.29 (6.40)
Median [Min, Max] 4.00 [0, 80.0]
Missing 1442 (3.1%)
Hours Exercise
Mean (SD) 5.78 (7.70)
Median [Min, Max] 4.00 [0, 80.0]
Missing 1442 (3.1%)
Hours Housework
Mean (SD) 10.3 (10.1)
Median [Min, Max] 8.00 [0, 168]
Missing 1442 (3.1%)
Hours Community
Mean (SD) 0.937 (2.87)
Median [Min, Max] 0 [0, 120]
Missing 1442 (3.1%)
Household Income
Mean (SD) 115000 (92100)
Median [Min, Max] 100000 [1.00, 3010000]
Missing 2940 (6.3%)
Male (binary)
0 29121 (62.8%)
1 17148 (37.0%)
Missing 108 (0.2%)
Neuroticism
Mean (SD) 3.49 (1.15)
Median [Min, Max] 3.50 [1.00, 7.00]
Missing 402 (0.9%)
Not Heterosexual Binary
0 42163 (90.9%)
1 3039 (6.6%)
Missing 1175 (2.5%)
NZ Deprevation Index 2018
Mean (SD) 4.77 (2.73)
Median [Min, Max] 4.00 [1.00, 10.0]
Missing 308 (0.7%)
NZSEI (Occupational Prestige Index)
Mean (SD) 54.1 (16.5)
Median [Min, Max] 54.0 [10.0, 90.0]
Missing 346 (0.7%)
Openness
Mean (SD) 4.96 (1.12)
Median [Min, Max] 5.00 [1.00, 7.00]
Missing 394 (0.8%)
Parent (binary)
0 13541 (29.2%)
1 32836 (70.8%)
Partner Binary
Mean (SD) 0.750 (0.433)
Median [Min, Max] 1.00 [0, 1.00]
Missing 535 (1.2%)
Political Conservative
Mean (SD) 3.59 (1.38)
Median [Min, Max] 4.00 [1.00, 7.00]
Missing 1902 (4.1%)
Religion Identification Level
Mean (SD) 2.38 (2.20)
Median [Min, Max] 1.00 [1.00, 7.00]
Missing 295 (0.6%)
Religion Church Binary
Mean (SD) 0.170 (0.376)
Median [Min, Max] 0 [0, 1.00]
Rural Gch 2018 Levels
High Urban Accessibility 28607 (61.7%)
Medium Urban Accessibility 8695 (18.7%)
Low Urban Accessibility 5639 (12.2%)
Remote 2581 (5.6%)
Very Remote 549 (1.2%)
Missing 306 (0.7%)
Right Wing Authoritarianism
Mean (SD) 3.28 (1.15)
Median [Min, Max] 3.20 [1.00, 7.00]
Missing 7 (0.0%)
Sample Frame Opt-In (binary)
0 45001 (97.0%)
1 1376 (3.0%)
Social Dominance Orientation
Mean (SD) 2.32 (0.961)
Median [Min, Max] 2.17 [1.00, 7.00]
Missing 1 (0.0%)
Smoker (binary)
0 41895 (90.3%)
1 3297 (7.1%)
Missing 1185 (2.6%)

Exposure Variable: Religious Service Attendance

Table 2 describes religious service attendance at baseline and during the exposure waves. Attendance was not collected in parts of wave 12 (2020–2021) and wave 13 (2021–2022). Because this exposure is central to our research question, we retained it and managed the missingness as follows: 1. Carry‑forward imputation. When a later wave recorded attendance, we filled the earlier gap with the most recent observed value and added an “imputed” flag. 2. Censoring. If no later attendance data existed, we treated the participant as censored from that point and applied inverse‑probability‑of‑censoring weights, as outlined in Method: Handling of missing data.

Including the exposure despite missingness preserves sample size, maintains comparability with previous NZAVS analyses, and avoids discarding information that could clarify causal pathways. The carry‑forward step may attenuate effects by introducing random measurement error, but it is unlikely to inflate them. Thus, any detected causal effects should be viewed as conservative estimates.

Table 2: Sample statistics for the exposure variable—religious service attendance across NZAVS time‑10 through time 14 cohort waves (2018 – 2023). Attendance was not recorded in 2020 and was only partially recorded in 2021.
Baseline Wave Exposure 1 Exposure 2 Exposure 3 Exposure 4
(N=46377) (N=38165) (N=31470) (N=31470) (N=27954)
Religious Attendance (weekly)
0 38479 (83.0%) 32263 (84.5%) 26812 (85.2%) 27136 (86.2%) 23850 (85.3%)
1 1517 (3.3%) 1038 (2.7%) 832 (2.6%) 698 (2.2%) 742 (2.7%)
2 1125 (2.4%) 847 (2.2%) 674 (2.1%) 636 (2.0%) 652 (2.3%)
3 907 (2.0%) 718 (1.9%) 583 (1.9%) 534 (1.7%) 543 (1.9%)
4 2478 (5.3%) 1916 (5.0%) 1490 (4.7%) 1425 (4.5%) 1319 (4.7%)
5 475 (1.0%) 352 (0.9%) 290 (0.9%) 271 (0.9%) 212 (0.8%)
6 326 (0.7%) 231 (0.6%) 170 (0.5%) 168 (0.5%) 139 (0.5%)
7 106 (0.2%) 79 (0.2%) 59 (0.2%) 52 (0.2%) 41 (0.1%)
8 964 (2.1%) 721 (1.9%) 560 (1.8%) 550 (1.7%) 456 (1.6%)

Outcome Variables

Table 3: Outcome variables measured at baseline (NZAVS time 10, years 2018-2019, and time 15, years 2023-2024).
Baseline Wave Outcome Wave
(N=46377) (N=46377)
Alcohol Frequency
Mean (SD) 2.16 (1.34) 2.04 (1.36)
Median [Min, Max] 2.00 [0, 5.00] 2.00 [0, 5.00]
Missing 1598 (3.4%) 25538 (55.1%)
Alcohol Intensity
Mean (SD) 2.17 (2.17) 2.00 (1.89)
Median [Min, Max] 2.00 [0, 30.0] 2.00 [0, 48.0]
Missing 2751 (5.9%) 27939 (60.2%)
Body Mass Index
Mean (SD) 27.2 (5.87) 27.8 (6.07)
Median [Min, Max] 26.2 [12.3, 73.6] 26.7 [13.2, 87.4]
Missing 1171 (2.5%) 25123 (54.2%)
Sleep
Mean (SD) 6.94 (1.13) 6.93 (1.12)
Median [Min, Max] 7.00 [2.50, 16.0] 7.00 [2.00, 16.0]
Missing 2365 (5.1%) 26025 (56.1%)
Hours of Exercise
Mean (SD) 5.78 (7.70) 6.42 (7.29)
Median [Min, Max] 4.00 [0, 80.0] 5.00 [0, 80.0]
Missing 1442 (3.1%) 25897 (55.8%)
Short Form Health
Mean (SD) 5.04 (1.17) 4.84 (1.17)
Median [Min, Max] 5.00 [1.00, 7.00] 5.00 [1.00, 7.00]
Missing 9 (0.0%) 25080 (54.1%)
Fatigue
Mean (SD) 1.64 (1.09) 1.62 (1.07)
Median [Min, Max] 2.00 [0, 4.00] 2.00 [0, 4.00]
Missing 500 (1.1%) 25117 (54.2%)
Anxiety
Mean (SD) 1.21 (0.773) 1.16 (0.764)
Median [Min, Max] 1.00 [0, 4.00] 1.00 [0, 4.00]
Missing 437 (0.9%) 25048 (54.0%)
Depression
Mean (SD) 0.585 (0.753) 0.526 (0.722)
Median [Min, Max] 0.333 [0, 4.00] 0.333 [0, 4.00]
Missing 440 (0.9%) 25048 (54.0%)
Rumination
Mean (SD) 0.847 (1.01) 0.759 (0.956)
Median [Min, Max] 1.00 [0, 4.00] 0 [0, 4.00]
Missing 538 (1.2%) 25108 (54.1%)
Body Satisfaction
Mean (SD) 4.24 (1.70) 4.23 (1.69)
Median [Min, Max] 4.00 [1.00, 7.00] 4.00 [1.00, 7.00]
Missing 537 (1.2%) 25931 (55.9%)
Forgiveness
Mean (SD) 5.02 (1.27) 5.20 (1.26)
Median [Min, Max] 5.33 [1.00, 7.00] 5.33 [1.00, 7.00]
Missing 11 (0.0%) 25125 (54.2%)
Perfectionism
Mean (SD) 3.18 (1.34) 2.98 (1.40)
Median [Min, Max] 3.00 [1.00, 7.00] 2.67 [1.00, 7.00]
Missing 6 (0.0%) 25174 (54.3%)
Self Control
Mean (SD) 4.40 (1.41) 4.48 (1.43)
Median [Min, Max] 4.50 [1.00, 7.00] 4.50 [1.00, 7.00]
Missing 23 (0.0%) 25310 (54.6%)
Self Esteem
Mean (SD) 5.14 (1.28) 5.27 (1.30)
Median [Min, Max] 5.33 [1.00, 7.00] 5.67 [1.00, 7.00]
Missing 400 (0.9%) 25104 (54.1%)
Sexual Satisfaction
Mean (SD) 4.58 (1.77) 4.38 (1.76)
Median [Min, Max] 5.00 [1.00, 7.00] 4.00 [1.00, 7.00]
Missing 2958 (6.4%) 26494 (57.1%)
Gratitude
Mean (SD) 5.90 (0.886) 5.90 (0.910)
Median [Min, Max] 6.00 [1.00, 7.00] 6.00 [1.00, 7.00]
Missing 9 (0.0%) 25102 (54.1%)
Life Satisfaction
Mean (SD) 5.29 (1.20) 5.22 (1.24)
Median [Min, Max] 5.50 [1.00, 7.00] 5.50 [1.00, 7.00]
Missing 223 (0.5%) 25348 (54.7%)
Meaning: Purpose
Mean (SD) 5.20 (1.42) 5.16 (1.41)
Median [Min, Max] 5.00 [1.00, 7.00] 5.00 [1.00, 7.00]
Missing 1224 (2.6%) 26377 (56.9%)
Meaning: Sense
Mean (SD) 5.71 (1.22) 5.78 (1.20)
Median [Min, Max] 6.00 [1.00, 7.00] 6.00 [1.00, 7.00]
Missing 160 (0.3%) 26039 (56.1%)
Personal Well-being Index
Mean (SD) 7.08 (1.66) 6.97 (1.74)
Median [Min, Max] 7.25 [0, 10.0] 7.25 [0, 10.0]
Missing 49 (0.1%) 25062 (54.0%)
Social Belonging
Mean (SD) 5.14 (1.08) 5.23 (1.10)
Median [Min, Max] 5.33 [1.00, 7.00] 5.33 [1.00, 7.00]
Missing 397 (0.9%) 25100 (54.1%)
Neighbourhood Belonging
Mean (SD) 4.19 (1.67) 4.63 (1.56)
Median [Min, Max] 4.00 [1.00, 7.00] 5.00 [1.00, 7.00]
Missing 252 (0.5%) 25906 (55.9%)
Social Support (perceived)
Mean (SD) 5.95 (1.12) 6.00 (1.13)
Median [Min, Max] 6.33 [1.00, 7.00] 6.33 [1.00, 7.00]
Missing 37 (0.1%) 25133 (54.2%)

S3: Confouding Control

Table 4: Table 4 presents single-world intervention graphs showing time-fixed and time-varying sources of bias in our six waves (baseline, four exposure waves, followed by the outcome wave.) Time-fixed confounders are included in the baseline wave. Time-varying confounders are included in each of the four treatment waves (abbreviated here by ‘\dots’ to declutter the graph). When there is more than one exposure wave, identifying causal effects requires adjustment for time-varying confounders (Robins and Hernan 2008; Bulbulia 2024; Richardson and Robins 2013).

For confounding control, we employ a modified disjunctive cause criterion (VanderWeele 2019), which involves:

  1. Identifying all common causes of both the treatment and outcomes.
  2. Excluding instrumental variables that affect the exposure but not the outcome.
  3. Including proxies for unmeasured confounders affecting both exposure and outcome.
  4. Controlling for baseline exposure and baseline outcome, serving as proxies for unmeasured common causes (VanderWeele, Mathur, and Chen 2020).

Additionally, we control for time-varying confounders at each exposure wave (Robins and Hernan 2008; Bulbulia 2024; Richardson and Robins 2013).

The covariates included for confounding control are described in Rosa et al. (2024).

Where there are multiple exposures, causal inference may be threatened by time-varying confounding (Bulbulia 2024).

S4: Causal Contrasts and Causal Assumptions

Notation

  • A_k: Observed religious service attendance at Wave k, for k = 1, \dots, 4.
  • Y_\tau: Outcome measured at the end of the study (Wave 5).
  • W_0: Confounders measured at baseline (Wave 0).
  • L_k: Time-varying confounders measured at Wave k (for k = 1, \dots, 4).

Shift Functions

Let \boldsymbol{\text{d}}(a_k)^+ represent the regular attendance treatment sequence and \boldsymbol{\text{d}}(a_k)^- the no attendance treatment sequence, where the interventions occure at each wave k = 1\dots 4; k\in \{0\dots 5\}. Formally:

Regular Attendance Regime \bigl(\boldsymbol{\text{d}}(a_k^+)\bigr)

\boldsymbol{\text{d}} (a_k^+) \;=\; \begin{cases} 4, & \text{if } A_k < 4,\\[6pt] A_k, & \text{otherwise.} \end{cases}

No Attendance Regime \bigl(\boldsymbol{\text{d}}(a_k^-)\bigr)

\boldsymbol{\text{d}}(a_k^-) \;=\; \begin{cases} 0, & \text{if } A_k > 0,\\[6pt] A_k, & \text{otherwise.} \end{cases}

Here, A_k is the observed attendance at Wave k. The shift function \boldsymbol{\text{d}} “nudges” A_k to a target level (four times per month or zero) only if the current value is below (for regular attendance) or above (for no attendance) that target. Across the four waves, these shifts form a sequence \boldsymbol{\bar{\boldsymbol{\text{d}}}}, which defines a complete intervention regime.

Causal Contrast

We compare the average well-being under the regular attendance regime, \boldsymbol{\bar{\boldsymbol{\text{d}}}}(a^+), to the average well-being under the no attendance regime, \boldsymbol{\bar{\boldsymbol{\text{d}}}}(a^-). Specifically, the average treatment effect (ATE) is given:

\text{ATE}^{\text{wellbeing}} \;=\; \mathbb{E} \Bigl[ Y_\tau\!\bigl(\boldsymbol{\text{d}}(a^+)\bigr) \;-\; Y_\tau\!\bigl(\boldsymbol{\text{d}}(a^-)\bigr) \Bigr].

Assumptions

To estimate this effect from observational data, we assume:

  1. Conditional Exchangeability: Once we condition on W_0 and each L_k, the interventions \boldsymbol{\bar{\boldsymbol{\text{d}}}}(a^+) or \boldsymbol{\bar{\boldsymbol{\text{d}}}}(a^-) are effectively random with respect to potential outcomes.
  2. Consistency: The potential outcome under a given regime matches the observed outcome when that regime is followed.
  3. Positivity: Everyone has a non-zero probability of receiving each level of attendance (i.e., a chance to be “shifted” up or down) given their covariates. The positivity assumption is the only causal assumption that can be evaluated with data. We evaluate this assumption in Appendix E).

Mathematically, for conditional exchangeability, we write:

\Bigl\{ Y\bigl(\boldsymbol{\text{d}}(a^+)\bigr), \; Y\bigl(\boldsymbol{\text{d}}(a^-)\bigr) \Bigr\} \coprod A_k | W_0, L_k That is, we assume the potential outcomes under each treatment regime are independent of each treatment at every time point, conditional on baseline confounders and time-varying confounders.

Under these assumptions, our statistical models permit us to estimate \text{ATE}^{\text{wellbeing}} from observational data. That contrast is (1) regularly attend religious service for four wave versus (2) never attend religious service at least four waves. We define the target population as the New Zealand Population from 2019-2024, the years in which measurements were taken.

S5: Year‑on‑year transitions in religious service attendance

These transition matrices capture shifts in states between consecutive waves. Each cell shows the count of individuals transitioning from one state to another. Rows are the initial state (From), columns the subsequent state (To). Diagonal entries (in bold) mark those who stayed in the same state.

Transition Matrix From Wave 2018 to Wave 2019
From / To State 0 State 1 State 2 State 3 State 4 State 5 State 6 State 7 State 8 Total
State 0 27321 405 174 71 126 26 13 8 68 28212
State 1 658 244 85 44 46 5 2 3 10 1097
State 2 239 105 199 104 96 12 13 2 21 791
State 3 114 54 110 171 173 18 8 4 15 667
State 4 155 71 127 205 897 124 64 16 91 1750
State 5 24 7 17 17 145 61 25 7 33 336
State 6 15 5 13 17 84 22 33 5 37 231
State 7 9 0 6 3 16 6 9 7 19 75
State 8 75 14 17 14 105 34 42 17 359 677

S6 Heterogeneity Treatment Effect Decision Flow

This following flowchart shows the decision logic:

    START: For each model
             |
             v
    STEP 1: EXCLUSION CHECK
    Is RATE QINI < 0 (stat sig) OR RATE AUTOC < 0 (stat sig)?
             |
        +----+----+
        |         |
       YES       NO
        |         |
        v         v
    EXCLUDED    STEP 2: RATE SELECTION CHECK
    (Stop)      Is RATE QINI > 0 (stat sig) OR RATE AUTOC > 0 (stat sig)?
                         |
                    +----+----+
                    |         |
                   YES       NO
                    |         |
                    v         v
                SELECTED   UNCLEAR
                    |         |
                    v         v
                STEP 3: QINI CURVE ANALYSIS
                (Applied to all non-excluded models)
                    |         |
                    v         v
                Already     Check: Is QINI curve positive
                selected     at any spend level?
                            |
                       +----+----+
                       |         |
                      YES       NO
                       |         |
                       v         v
                   SELECTED   UNCLEAR
                             (final)

References

Altemeyer, Bob. 1996. The Authoritarian Spectre. London: Harvard University Press.
Atkinson, J, C Salmond, and P Crampton. 2019. “NZDep2018 Index of Deprivation, User’s Manual.” Wellington.
Bulbulia, J. A. 2024. “Methods in Causal Inference Part 2: Interaction, Mediation, and Time-Varying Treatments.” Evolutionary Human Sciences 6: e41. https://doi.org/10.1017/ehs.2024.32.
Bulbulia, J. A., J. H. Shaver, L. Greaves, R. Sosis, and C. G. Sibley. 2015. “Religion and Parental Cooperation: An Empirical Test of Slone’s Sexual Signaling Model.” In The Attraction of Religion: A Sexual Selectionist Account, edited by Slone J. D. amd Van Slyke J., 29–62. Bloomsbury Press.
Darwin, Charles. 1872. The Descent of Man, and Selection in Relation to Sex. Vol. 2. D. Appleton.
Dew, Jeffrey P, Jeremy E Uecker, and Brian J Willoughby. 2020. “Joint Religiosity and Married Couples’ Sexual Satisfaction.” Psychology of Religion and Spirituality 12 (2): 201.
Fahy, Katie M., Alan Lee, and Barry J. Milne. 2017. New Zealand Socio-Economic Index 2013. Wellington, New Zealand: Statistics New Zealand-Tatauranga Aotearoa.
Fraser, Gloria, Joseph Bulbulia, Lara M. Greaves, Marc S. Wilson, and Chris G. Sibley. 2020. “Coding Responses to an Open-Ended Gender Measure in a New Zealand National Sample.” The Journal of Sex Research 57 (8): 979–86. https://doi.org/10.1080/00224499.2019.1687640.
Greaves, Lara M, Fiona Kate Barlow, Carol HJ Lee, Correna M Matika, Weiyu Wang, Cinnamon-Jo Lindsay, Claudia JB Case, et al. 2017. “The Diversity and Prevalence of Sexual Orientation Self-Labels in a New Zealand National Sample.” Archives of Sexual Behavior 46: 1325–36.
Hagerty, Bonnie M. K., and Kathleen Patusky. 1995. “Developing a Measure Of Sense of Belonging:” Nursing Research 44 (1): 9–13. https://doi.org/10.1097/00006199-199501000-00003.
Health, Ministry of. 2013. “The New Zealand Health Survey: Content Guide 2012-2013.” Princeton University Press.
McCullough, Michael E, Shelley D Kilpatrick, Robert A Emmons, and David B Larson. 2001. “Is Gratitude a Moral Affect?” Psychological Bulletin 127 (2): 249.
Richardson, Thomas S., and James M. Robins. 2013. “Single World Intervention Graphs: A Primer.” In. Citeseer. https://core.ac.uk/display/102673558.
Robins, James, and Miguel Hernan. 2008. “Estimation of the Causal Effects of Time-Varying Exposures.” Chapman & Hall/CRC Handbooks of Modern Statistical Methods, 553–99.
Rosa, Pedro A de la, Richard G Cowden, Joseph A Bulbulia, Chris G Sibley, and Tyler J VanderWeele. 2024. “Effects of Screen-Based Leisure Time on 24 Subsequent Health and Wellbeing Outcomes: A Longitudinal Outcome-Wide Analysis.” International Journal of Behavioral Medicine, 1–20. https://doi.org/https://doi.org/10.1007/s12529-024-10307-0.
Shaver, John H, Radim Chvaja, Laure Spake, Anushé Hassan, Jainaba Badjie, Andrew M Prentice, Carla Cerami, Rebecca Sear, Mary K Shenk, and Richard Sosis. 2024. “Religious Involvement Is Associated with Higher Fertility and Lower Maternal Investment, but More Alloparental Support Among Gambian Mothers.” American Journal of Human Biology, e24144.
Shaver, John H, Chris G Sibley, Richard Sosis, Deane Galbraith, and Joseph Bulbulia. 2019. “Alloparenting and Religious Fertility: A Test of the Religious Alloparenting Hypothesis.” Evolution and Human Behavior 40 (3): 315–24.
Sibley, Chris G. 2021. “Sampling Procedure and Sample Details for the New Zealand Attitudes and Values Study.” https://doi.org/10.31234/osf.io/wgqvy.
Sibley, Chris G, Nils Luyten, Missy Purnomo, Annelise Mobberley, Liz W Wootton, Matthew D Hammond, Nikhil Sengupta, et al. 2011. “The Mini-IPIP6: Validation and Extension of a Short Measure of the Big-Six Factors of Personality in New Zealand.” New Zealand Journal of Psychology 40 (3): 142–59.
Sidanius, Jim, and Felicia Pratto. 1999. Social Dominance: An Intergroup Theory of Social Hierarchy and Oppression. Cambridge: Cambridge University Press.
Statistics New Zealand. 2017. Statistical Standard for Geographic Areas 2018 (SSGA18). Wellington, New Zealand: Statistics New Zealand. https://www.stats.govt.nz/methods/statistical-standard-for-geographic-areas-2018/.
VanderWeele, Tyler J. 2019. “Principles of Confounder Selection.” European Journal of Epidemiology 34 (3): 211–19.
VanderWeele, Tyler J, Maya B Mathur, and Ying Chen. 2020. “Outcome-Wide Longitudinal Designs for Causal Inference: A New Template for Empirical Studies.” Statistical Science 35 (3): 437–66.
Verbrugge, Lois M. 1997. “A Global Disability Indicator.” Journal of Aging Studies 11 (4): 337–62. https://doi.org/10.1016/S0890-4065(97)90026-8.
Whitehead, Jesse, Gabrielle Davie, Brandon de Graaf, Sue Crengle, Ross Lawrenson, Rory Miller, and Garry Nixon. 2023. “Unmasking Hidden Disparities: A Comparative Observational Study Examining the Impact of Different Rurality Classifications for Health Research in Aotearoa New Zealand.” BMJ Open 13 (4): e067927.